load('BloombergData20240205Govt.mat')
res=zeros(days1+days2+1,29,2);
for ev=1:29
    if ev<=8
        bo=302:320; %Danish bonds
    else
        bo=1:74; %German bonds
    end
    %If all bonds have a yield of zero on a given day, it is because the
    %observations are missing (i.e. non-trading day)
    for i=1:length(AllDates)
        ind(i)=sum(abs(bonddat(i,2,bo)))>0;
    end
    ind=find(ind);
    AllDates2=AllDates(ind);
    bonddat2=bonddat(ind,:,:);
    bo2=101:142; %Dutch control bonds
    if ev==1
        ebond=305;%20-Jan-2022 first day of Danish GB trading. 
        eday=min(find(bonddat2(:,2,ebond)~=0));
    elseif ev==2
        ebond=305;
        eday=find(AllDates2==datenum('17March2022'));
    elseif ev==3
        ebond=305;
        eday=find(AllDates2==datenum('19May2022'));
    elseif ev==4
        ebond=305;
        eday=find(AllDates2==datenum('07July2022'));
    elseif ev==5
        ebond=305;
        eday=find(AllDates2==datenum('08Sept2022'));
    elseif ev==6
        ebond=305;
        eday=find(AllDates2==datenum('24Oct2022'));
    elseif ev==7
        ebond=316;
        eday=find(AllDates2==datenum('27Sept2023'));
    elseif ev==8
        ebond=316;
        eday=find(AllDates2==datenum('09Nov2023'));
    elseif ev==9
        ebond=62;
        eday=min(find(bonddat2(:,2,ebond)~=0));
    elseif ev==10
        ebond=62;
        eday=find(AllDates2==datenum('25Jan2023'));
    elseif ev==11
        ebond=62;
        eday=find(AllDates2==datenum('08June2023'));
    elseif ev==12
        ebond=22;%German 3
        eday=min(find(bonddat2(:,2,ebond)~=0)); 
    elseif ev==13
        ebond=22;
        eday=find(AllDates2==datenum('03March2022'));
    elseif ev==14
        ebond=22;
        eday=find(AllDates2==datenum('21July2022'));
    elseif ev==15
        ebond=18;
        eday=min(find(bonddat2(:,2,ebond)~=0)); 
    elseif ev==16
        ebond=18;
        eday=find(AllDates2==datenum('21October2021'));
    elseif ev==17
        ebond=18;
        eday=find(AllDates2==datenum('05May2022'));
    elseif ev==18
        ebond=18;
        eday=find(AllDates2==datenum('02Nov2022'));
   elseif ev==19
        ebond=3;%German 1
        eday=min(find(bonddat2(:,2,ebond)~=0)); 
    elseif ev==20
        ebond=3;
        eday=find(AllDates2==datenum('02June2022'));   
    elseif ev==21
        ebond=73;
        eday=min(find(bonddat2(:,2,ebond)~=0)); 
    elseif ev==22
        ebond=73;
        eday=find(AllDates2==datenum('22March2023'));
     elseif ev==23
        ebond=73;
        eday=find(AllDates2==datenum('31Aug2023')); 
     elseif ev==24
        ebond=73;
        eday=find(AllDates2==datenum('24January2024')); 
    elseif ev==25
        ebond=68;
        eday=min(find(bonddat2(:,2,ebond)~=0)); 
    elseif ev==26
        ebond=68;
        eday=find(AllDates2==datenum('1Nov2023'));
    elseif ev==27
        ebond=68;
        eday=find(AllDates2==datenum('24Jan2024'));
   elseif ev==28
        ebond=69;
        eday=min(find(bonddat2(:,2,ebond)~=0)); 
    elseif ev==29
        ebond=69;
        eday=find(AllDates2==datenum('06July2023'));   
    else end

    [C,ia,ic] = unique(bondinfo(bo2,4)); %4=bond maturity
    bo2=sort(bo2(ia)); 
    mat=(bondinfo(bo,4)-AllDates2(eday))/365.25; %OBS mat has descending maturities
    mat2=(bondinfo(bo2,4)-AllDates2(eday))/365.25; %OBS mat2 has descending maturities
    minmat=1/12;
    bo=bo(mat>minmat&squeeze(bonddat2(eday,yno,bo)~=0));   
    outamount=bondinfo(bo,1); %1=outstanding amount
    mat=(bondinfo(bo,4)-AllDates2(eday))/365.25;
    bo2=bo2(mat2>minmat&squeeze(bonddat2(eday,yno,bo2)~=0)); 
    mat2=(bondinfo(bo2,4)-AllDates2(eday))/365.25;
    bonddat3=bonddat2./0; %bonddat3 has spread relative to control bonds
    for i=1:length(AllDates2)
        for b=1:length(bo)
            bo1=min(find(mat2<mat(b)));
            [M I]=min(abs(mat2-mat(b)));
            bonddat3(i,yno,bo(b))=bonddat2(i,yno,bo(b))-bonddat2(i,yno,bo2(I));
        end
    end
    av=zeros(length(AllDates2),2)./0;
    for i=1:length(AllDates2)
        av(i,1)=sum(squeeze(bonddat3(i,yno,bo)).*outamount)/sum(outamount);
        indblack=ismember(bo,[305 62 22 18 3])==0;
        av(i,2)=sum(squeeze(bonddat3(i,yno,bo(indblack))).*outamount(indblack))/sum(outamount(indblack));
    end 
    res(:,ev,:)=100*av(eday-days1:eday+days2,:);
end


